function [u,varargout]=UtilityFlow_Sep(NU_C,NU_L,PSI_L,zw,input,InputType)
%% Preliminaries
if isscalar(zw) && isvector(input)
    zw      =   zw*ones(size(input));
end

if nargin<=5
    InputType   =   'Consumption';
end

%% Solve the Optimal Labor Supply
switch InputType
    case 'Consumption'
        c           =   input;
        % Labor Supply
        if NU_C<0 || NU_L<=0 || PSI_L<=0
            error('Undertermined Labor Supply');
        end
        l           =   ( zw.*(c.^(-NU_C))/PSI_L ).^( 1/NU_L );
        chat        =   c-zw.*l;
        varargout{1}=   chat;
        varargout{2}=   l;
    case 'Expenditure'
        chat        =   input;
        if NU_C<=0 || NU_L<0 || PSI_L<=0
            error('Undertermined Labor Supply');
        elseif NU_L==0
            c           =   (zw/PSI_L).^(1/NU_C);
        else
            Temp        =   ( zw.^(1+1/NU_L) )*( PSI_L^(-1/NU_L) );
            TempInd     =   Temp>1e-6;
            % Corresponding to the Unemployment
            c           =   chat;
            % Corresponding to the Normal Employed State
            if any(TempInd)
                temp        =   Temp(TempInd);
                temp_chat   =  chat(TempInd);
                TempFun     =   @(c) c - temp.*( c.^(-NU_C/NU_L) ) - temp_chat;
                c0          =   temp.^( NU_L/(NU_C+NU_L) );
                lb          =   max(temp_chat,0);
                ub          =   lb + 2*c0;
                c(TempInd)  =   VecFun_Bisection(TempFun,lb,ub,1e-7);
            end
        end
        l           =   (c-chat)./zw;
        varargout{1}=   c;
        varargout{2}=   l;
end

%% Utility
if NU_C==1
    u       =   log(c) - PSI_L*( l.^(1+NU_L)/(1+NU_L) );
else
    u       =   c.^(1-NU_C)/(1-NU_C) ...
                - PSI_L*( l.^(1+NU_L)/(1+NU_L) );
end
